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Abstract 

We describe an exact approach for calculating transition probabilities 
and waiting times in finite-state discrete-time Markov processes. All the 
states and the rules for transitions between them must be known in ad- 
vance. We can then calculate averages over a given ensemble of paths for 
both additive and multiplicative properties in a non-stochastic and non- 
iterative fashion. In particular, we can calculate the mean first-passage 
time between arbitrary groups of stationary points for discrete path sam- 
pling databases, and hence extract phenomenological rate constants. We 
present a number of examples to demonstrate the efficiency and robustness 
of this approach. 



*E-mail: sat39@caiii.ac.uk 
^E-mail: dw34@cam.ac.uk 



1 



1 Introduction 



Stochastic processes are widely used to treat phenomena with random factors and 
noise. Markov processes are an important class of stochastic processes for which 
future transitions do not depend upon how the current state was reached. Markov 
processes restricted to a discrete, finite, or countably infinite state space are called 
Markov chains. Many interesting problems of chemical kinetics concern the 
analysis of finite-state samples of otherwise infinite state space. ^ 

When analysing the kinetic databases obtained from discrete path sampling 
(DPS) studies^ it can be difficult to extract the phenomenological rate constants 
for processes that occur over very long time scales.^ DPS databases are composed 
of local minima of the potential energy surface (PES) and the transition states 
that connect them. While minima correspond to mechanically stable structures, 
the transition states specify how these structures interconvert and can be used to 
calculate the corresponding rates. Whenever the potential energy barrier for the 
event of interest is large in comparison with ksT the event becomes rare, where 
T is the temperature and is Boltzmann's constant. 

The most important tools previously employed to extract kinetic informa- 
tion from a DPS stationary point database are the master equation,^ kinetic 
Monte Carlo^'^ (KMC) and matrix multiplication (MM) methods.^ The system 
of linear master equations in its matrix formulation can be solved numerically 
to yield the time evolution of the occupation probabilities starting from an arbi- 
trary initial distribution. This approach works well only for small problems, as 
the diagonalisation of the transition matrix, P, scales as the cube of the num- 
ber of states.^ In addition, numerical problems arise when the magnitude of the 
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eigenvalues corresponding to the slowest relaxation mode approaches the preci- 
sion of the zero eigenvalue corresponding to equilibrium.^ The KMC approach is 
a stochastic technique that is commonly used to simulate the dynamics of various 
physical and chemical systems, examples being formation of crystal structures/" 
nanoparticle growth^^ and diffusion. The MM approach provides a way to sum 
contributions to phenomenological two-state rate constants from pathways that 
contain progressively more steps. It is based upon a steady-state approximation, 
and provides the corresponding solution to the linear master equation.^' The 
MM approach has been used to analyse DPS databases in a number of systems 
ranging from Lennard- Jones clusters^' to biomolecules.^^'^^ 

Both the standard KMC and MM formulations provide rates at a computa- 
tional cost that generally grows exponentially as the temperature is decreased. 
In this contribution we describe alternative methods that are deterministic and 
formally exact, where the computational requirements are independent of the 
temperature and the time scale on which the process of interest takes place. 

1.1 Graph Theory Representation of a Finite-state Markov 
Chain 

To fully define a Markov chain it is necessary to specify all the possible states of 
the system and the rules for transitions between them. Graph theoretical repre- 
sentations of finite-state Markov chains are widely used.^'^''"^^ Here we adopt a 
digraph^"'^^ representation, where nodes represent the states and edges represent 
the transitions with non-zero probabilities. The edge ejj describes a transition 
from node j to node i and has a probability Pjj associated with it, which is com- 
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monly known as a routing or branching probability. A node can be connected to 
any number of other nodes. Two nodes of a graph are adjacent if there is an edge 
between them.^^ 

For digraphs all connections of a node are classified as incoming or outgoing. 
The total number of incoming connections is the in-degree of a node, while the 
total number of outgoing connections is the out-degree. In a symmetric digraph 
the in-degree and out-degree are the same for every node.^^ Adjln[i] is the set of 
indices of all nodes that are connected to node i via incoming edges that finish 
at node i. Similarly, AdjOut[i] is the set of indices of all the nodes that are 
connected to node i via outgoing edges from node i. The degree of a graph is the 
maximum degree of all of its nodes. The expectation value for the degree of an 
undirected graph is the average number of connections per node. 

For any node i the transition probabilities Pj^i add up to unity. 



where the sum is over all j G AdjOut[i]. Unless specified otherwise all sums are 
taken over the set of indices of adjacent nodes or, since the branching probability 
is zero for non-adjacent nodes, over the set of all nodes. 

In a computer program dense graphs are usually stored in the form of adja- 
cency matrices. For sparse graphs^° a more compact but less efficient adjacency- 
lists-based data structure exists. To store a graph representation of a Markov 
chain, in addition to connectivity information (available from the adjacency ma- 
trix), the branching probabilities must be stored. Hence for dense graphs the 
most convenient approach is to store a transition probability matrix^ with transi- 
tion probabilities for non-existent edges set to zero. For sparse graphs, both the 
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adjacency list and a list of corresponding branching probabilities must be stored. 
1.2 The Kinetic Monte Carlo Method 

The KMC method can be used to generate a memoryless (Markovian) random 
walk and hence a set of trajectories connecting initial and final states in a DPS 
database. Many trajectories are necessary to collect appropriate statistics. Ex- 
amples of pathway averages that are usually obtained with KMC are the mean 
path length and the mean first-passage time. Here the KMC trajectory length 
is the number of states (local minima of the PES in the current context) that 
the walker encounters before reaching the final state. The first-passage time is 
defined as the time that elapses before the walker reaches the final state. For a 
given KMC trajectory the first-passage time is calculated as the sum of the mean 
waiting times in each of the states encountered. 

An efficient way to propagate KMC trajectories was suggested by Bortz, Ka- 
los, and Lebowitz (BKL).^ According to the BKL algorithm, a step is chosen 
in such a way that the ratios between transition probabilities of different events 
are preserved, but rejections are eliminated. Fig. 1 explains this approach for a 
simple discrete-time Markov chain. The evolution of an ordinary KMC trajectory 
is monitored by the 'time' parameter n G W, which actually corresponds to the 
number of steps. ^ The random walker is in state 1 at 'time' n = 0. The KMC 
trajectory is terminated whenever an absorbing state is encountered. As Pi^i 
approaches unity transitions out of state 1 become rare. To ensure that every 
time a random number is generated (one of the most time consuming steps in a 
KMC calculation) a move is made to a neighbouring state we average over the 
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transitions from state 1 to itself to obtain the Markov chain depicted in Fig. 1 (b). 
Transitions from state 1 to itself can be modelled by a Bernoulli process^^ with 
the probability of success equal to Pi^i. The average 'time' for escape from state 
1 is obtained as 

oo ^ 

n = (1 - + i)(Pi,ir = nzrjr-y (2) 

n=0 ^'^ 

which can be used as a measure of the efficiency of trapping. In the BKL 
scheme, transition probabilities out of state 1 are renormalised as: 

Similar ideas underlie the accelerated Monte Carlo algorithm suggested by Novotny.^^ 
Both the BKL and Novotny's methods can be many orders of magnitude faster 
than the standard KMC method when kinetic traps are present. 

In chemical kinetics transitions out of a state are usually described using a 
Poisson process, which can be considered a continuous-time analogue of Bernoulli 
trials. The transition probabilities are determined from the rates of the underlying 
transitions as 

= v^' (4) 



where ki j is the rate constant for transitions from j to i, etc. There may be 
several escape routes from a given state. Transitions from any state to directly 
connected states are treated as competing independent Poisson processes, which 
together generate a new Poisson distribution.^^ n independent Poisson processes 
with rates ki, /c2, /ca, . . . , kn combine to produce a Poisson process with rate 
k = Yl^=i ^i- The waiting time for a transition to occur to any connected state 
is then exponentially distributed as /c exp(— /ct).^^ 
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Given the exponential distribution of waiting times the mean waiting time 



is simply rf. When the average of the distribution of times is the property of 
interest, and not the distribution itself, it is sufficient to increment the simulation 
time by the mean waiting time rather than by a value drawn from the appropriate 
distribution.^' This modification to the original KMC formulation^^' reduces 
the cost of the method and accelerates the convergence of KMC averages without 
affecting the results. 

1.3 Discrete Path Sampling 

The result of a DPS simulation is a database of local minima and transition 
states from the PES.^'^'^^ To extract thermodynamic and kinetic properties from 
this database we require partition functions for the individual minima and rate 
constants, /cq,,/3, for the elementary transitions between adjacent minima (3 and 
a. We usually employ harmonic densities of states and statistical rate theory 
to obtain these quantities, but these details are not important here. To analyse 
the global kinetics we further assume Markovian transitions between adjacent 
local minima, which produces a set of linear (master) equations that governs the 
evolution of the occupation probabilities towards equilibrium^' 



where Pa{t) is the occupation probability of minimum a at time t. 

All the minima are classified into sets A, B and /, where A and B are the 
two states of interest and / corresponds to 'intervening' minima. When local 



in state i before escape, Tj, is ^/Ylj^jA^ variance of the waiting time 
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equilibrium is assumed within the A and B sets we can write 

Pa{t) = ^^^^ and n(t) = S^, (6) 

where -Pa(^) = SaG^-^a(^) and PBif) = J^beB^bit)- If the steady-state approxi- 
mation is apphed to all the intervening states z G /, so that 

dPi{t) 



dt 



0, (7) 



then Eq. (5) can be written as^ 

dPA{t) 



kA,BPB{t) - kB,APA{t), 



dt 

= kB,APA{t)-kA,BPB{t). 

The rate constants kA,B and kB,A for forward and backward transitions between 
states A and B are the sums over all possible paths within the set of interven- 
ing minima of the products of the branching probabilities corresponding to the 
elementary transitions for each path: 

h ■ h - ■ h - ■ h- , P^^ 

'^A,B - 2^ 



^\^k ■ ■ S^k ■ P^R 



(9) 



peq ) 



^ ] Pa,iiPii,i2 ' ' ' Pin-i,i 



B 

and similarly for ks^A-^ The superscript 'SS' specifies that the DPS rate constant 
formula was derived employing the steady-state approximation, as in the previous 
versions of the DPS method.^' The sum is over all paths that begin from a state 
6 G -B and end at a state a G A, and the prime indicates that paths are not allowed 
to revisit states in B. In previous contributions^' ^^"^^ this sum was evaluated 
using a weighted adjacency matrix multiplication method. The contributions 
of individual discrete paths to the total rate constants were also calculated in 



this way. However, analytic results from the theory of one-dimensional random 
walk^^"^'^ are now employed instead. It is also possible to evaluate rates without 
invoking the steady-state approximation,^^ as discussed in the following sections. 

1.4 KMC and Steady-state Averages 

We now show that the evaluation of the steady-state sum in Eq. (9) and the 
calculation of KMC averages are two closely related problems. 

For KMC simulations we define sources and sinks that coincide with the set of 
initial states, B, and final states. A, respectively. Every cycle of KMC simulation 
involves the generation of a single KMC trajectory connecting a node b ^ B and 
a node a & A. A source node b is chosen from set B with probability P^'^/P^. 

We can formulate the calculation of the mean first-passage time from B to A 
in graph theoretical terms as follows. Let the digraph consisting of nodes for all 
local minima and edges for each transition state be Q. The digraph consisting of 
all nodes except those belonging to region A is denoted by G. We assume that 
there are no isolated nodes in Q, so that all the nodes in A can be reached from 
every node in G, in one or more steps. Suppose we start a KMC simulation from 
a particular node (3 E G. Let Pa{n) be the expected occupation probability of 
node a after n KMC steps, with initial conditions Ppi^O) = 1 and Pa^fsiO) = 0. 
We further define an escape probability for each a G G as the sum of branching 
probabilities to nodes in A, i.e. 



KMC trajectories terminate when they arrive at an A minimum, and the expected 
probability transfer to the A region at the nth KMC step is '^aeG^a^a{n)- 




(10) 
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If there is at least one escape route from G to A with a non-zero branching 
probabihty, then eventually all the occupation probabilities in G must tend to 
zero and 

oo 

n=0 aGG 

for any (3 E G. We now rewrite Pa{n) as a sum over all n-step paths that start 
from j3 and end at a without leaving G. Each path contributes to Pain) according 
to the appropriate product of n branching probabilities, so that 

oo 

o6G n=0 
oo 

~ ^ ^ ^ ^ ^ ^ Pa,in-iPin-.i,i„-2 " " " -^i2,i\-^h,P (12) 
aeG n=0 S(n) 

aeG 

where 'E{n) denotes the set of n-step paths that start from (3 and end at a without 
leaving G, and the last line defines the pathway sum S'^^. 

It is clear from the last line of Eq. (12) that for fixed /? the quantities S^S^^ 
define a probability distribution. However, the pathway sums are not proba- 
bilities, and may be greater than unity. In particular, S^^ ^ 1 because the path 
of zero length is included, which contributes one to the sum. Furthermore, the 
normalisation condition on the last line of Eq. (12) places no restriction on S^^ 
terms for which vanishes. 

We can also define a probability distribution for individual pathways. Let 
be the product of branching probabilities associated with a path C, so that 

oo 

Sa,^ = E E = E (13) 

where a ^ j3 is the set of all appropriate paths from /5 to a of any length that 
can visit and revisit any node in G. If we focus upon paths starting from minima 
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in region B 

beB ^ aeG £,€a^b beB ^ o^Ga ^€a^b 

where is the set of nodes in G that are adjacent to A minima in the complete 
graph Q, since vanishes for all other nodes. We can rewrite this sum as 

peq pcq 

E 4f£fw,= E ^>^«= E n = i. (15) 

which defines the non-zero pathway probabilities for all paths that start from 
any node in set B and finish at any node in set A. The paths C, ^ A ^ B can 
revisit any minima in the G set, but include just one A minimum at the terminus. 
Note that and can be used interchangeably if there is only one state in set 
B. 

The average of some property, Q^, defined for each KMC trajectory, ^, can be 
calculated from the as 

^eA<~B 

Of course, KMC simulations avoid this complete enumeration by generating tra- 
jectories with probabilities proportional to V^, so that a simple running average 
can be used to calculate (Q^)- In the following sections we will develop alternative 
approaches based upon evaluating the complete sum, which become increasingly 
efficient at low temperature. We emphasise that these methods are only applica- 
ble to problems with a finite number of states, which are assumed to be known 
in advance. 

The evaluation of the sum based on the steady-state approximation defined 
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in Eq. (9) can also be rewritten in terms of pathway probabilities: 



oo 




,eq 
B 



n=0 ■E{n) 



oo 



peq 




n=0 S(n) 



(17) 



where the prime on the summation indicates that the paths are not allowed to 
revisit the B region. We have also used the fact that ki^^i, = Pi„^b/Tb- 

A digraph representation of the restricted set of pathways defined in Eq. (17) 
can be created if we allow sets of sources and sinks to overlap. In that case all the 
nodes AU B are defined to be sinks and all the nodes in B are sources, i.e. every 
node in set B is both a source and a sink. The required sum then includes all 
the pathways that finish at sinks of type A, but not those that finish at sinks 
of type B. This situation, where sets of sources and sinks (partially) overlap, is 
discussed in detail in §4. 

1.5 Mean Escape Times 

Since the mean first-passage time between states B and A, or the escape time 
from a subgraph, is of particular interest, we first illustrate a means to derive 
formulae for these quantities in terms of pathway probabilities. 

The average time taken to traverse a path ^ = ai, 0:2, as, . . . , ai(^) is calculated 
as = Tq,^ + Tq,2 + r^g, . . . , rQ,j^^j_j, where is the mean waiting time for escape 
from node a, as above, ak identifies the kth node along path ^, and 1{C,) is the 
length of path ^. The mean escape time from a graph G if started from node (3 
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is then 



E (18) 

If we multiply every branching probability, -Po,/3, that appears in by exp(^ra) 
then the mean escape time can be obtained as: 



n-G 



d 
dc 

d_ 



"i(€)'"i(€)-i ai(e)_i,ai(5)_2'= 



P e 



C=o 



E -PQi(«)'"i(«)-i-^°i(€)-i'^ 



■i(4)-2 • • • 02,01'-' 



C=o 



(19) 



This approach is useful if we have analytic results for the total probability S^, 
which may then be manipulated into formulae for T^^, and is standard practice 
in probability theory.^^'^^ The quantity Pa,f3^'^'^'^ is similar to the X probabil- 
ity' described in Ref. 35. Analogous techniques are usually employed to obtain 
and higher moments of the first-passage time distribution from analytic ex- 
pressions for the first-passage probability generating function (see, for example, 
Refs. 31,32). We now define Po.,p = Pa,^^''^'^ and the related quantities 



S. 



G 



p = cG (T^ 
/ , a,o <~-OL ^ 1 



'^o,/3 



P P 



and 



a£G 



(20) 
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Note that 



G 



etc., while the mean escape time can now be written as 



C=o 



dc 



(21) 



J c=o 



In the remaining sections we show how to calculate the pathway probabilities, 
V^, exactly, along with pathway averages, such as the waiting time. 



2 Complete graphs 

In a complete digraph each pair of nodes is connected by two oppositely directed 
edges. The complete graph with graph nodes is denoted K^, and has 
nodes and A^(A^ — 1) edges, remembering that we have two edges per connection. 
Due to the complete connectivity we need only consider two cases: when the 
starting and finishing nodes are the same and when they are distinct. We employ 
complete graphs for the purposes of generality. An arbitrary graph G^v is a 
subgraph of Kn with transition probabilities for non-existent edges set to zero. 
All the results in this section are therefore equally applicable to arbitrary graphs. 

a,l3 



The Sj^g can be derived analytically: 



(22) 



^1,1 = Yl + Wi,3,l + Wi,2,3,l + Wi,3,2,l) Yl ("^2,3,2)' 

n=0 \ m=0 

1 - >V2,3,2 

1 _ >V^ - >V2,3,2 - >Vi,3,l - Wi,2,3,l - Wi,3,2,l ' 

00 

S2J = XI ("^2,3,2)" (P2,l+W2,3,l)5ff 
n=0 

A,l + >V2,3,1 

1 _ 2,1 - >V2,3,2 - >Vi,3,l - >Vi,2,3,l - >Vi,3,2,l ' 

where, as before, VVi^2,i = -^1,2-^2,1; etc. The results for any other possibility can 
be obtained by permuting the node indices appropriately. 
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Pathway sums for larger complete graphs can be obtained by recursion. For 
S^lf any path leaving from and returning to can be broken down into a step 
out of to any i < N, all possible paths between i and j < N — 1 within Kn^i, 
and finally a step back to from j. All such paths can be combined together in 
any order, so we have a multinomial distribution:^® 

oo /N~l /N~l 

n=0 \i=l \j=l 

(N-l N-l 
i=l j=l 

To evaluate S^^ we break down the sum into all paths that depart from and 
return to A^, followed by all paths that leave node A^ and reach node 1 without 
returning to A^. The first contribution corresponds to a factor of S^'^, and the 
second produces a factor Pi^]\fS^f~^ for every i < N: 

N-l 
i=l 

where S^i is defined to be unity. Any other S^'^ can be obtained by a permuta- 
tion of node labels. 

Algorithm 1 provides an implementation of the above formulae optimised for 
incomplete graphs. The running time of Algorithm 1 depends strongly on the 
graph density. (A digraph in which the number of edges is close to the maximum 
value of A^(A^ — 1) is termed a dense digraph. For the algorithm runs 
in 0{N'^^) time, while for an arbitrary graph it scales as 0{d?^), where d is 
the average degree of the nodes. For chain graphs the algorithm runs in 0{N) 
time and has linear memory requirements. For complete graphs an alternative 
implementation with (9((A^!)^) scaling is also possible. 
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Although the scahng of the above algorithm with may appear disastrous, 
it does in fact run faster than standard KMC and MM approaches for graphs 
where the escape probabilities are several orders of magnitude smaller than the 
transition probabilities (Algorithm 1). Otherwise, for anything but moderately 
branched chain graphs. Algorithm 1 is significantly more expensive. However, the 
graph-transformation-based method presented in §2.1 yields both the pathway 
sums and the mean escape times for a complete graph Kj^ in 0{N^) time, and 
is the fastest approach that we have found. 

Mean escape times for are readily obtained from the results in equation 
(22) by applying the method outlined in §1.5: 

^K, ^ ri(l - W2,3,2) + r2(P2,l + W2,3,l) + T^jP^^i + ^3,2,1) ^25) 

We have verified this result numerically for various values of the parameters 
and Pci^0 and obtained quantitative agreement. Fig. 2 demonstrates how the 
advantage of exact summation over KMC and MM becomes more pronounced as 
the escape probabilities become smaller. 

2.1 Mean escape time from 

The problem of calculating the properties of a random walk on irregular net- 
works was addressed previously by Goldhirsch and Gefen."^^'^^ They described a 
generating-function-based method where an ensemble of pathways is partitioned 
into 'basic walks'. To the best of our knowledge only one"^^ out of the 30 pa- 
pgj,g3i,32,36,39-65 ^]-^g work of Goldhirsch and Gefen^^ is an application, 

perhaps due to the complexity of the method. Here we present a graph trans- 
formation (GT) approach for calculating the pathway sums and the mean escape 
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times for K^^. In general, it is applicable to arbitrary digraphs, but the best 
performance is achieved when the graph in question is dense. A sparse-optimised 
version of the GT method will be discussed in §3. 

The GT approach is similar in spirit to the ideas that lie behind the mean value 
analysis and aggregation/disaggregation techniques commonly used in the per- 
formance and reliability evaluation of queueing networks. ^'^^"^^ It is also loosely 
related to dynamic graph algorithms, ^^"''^ which are used when a property is cal- 
culated on a graph subject to dynamic changes, such as deletions and insertions 
of nodes and edges. The main idea is to progressively remove nodes from a graph 
whilst leaving the average properties of interest unchanged. For example, sup- 
pose we wish to remove node x from graph G to obtain a new graph G' . Here we 
assume that x is neither source nor sink. Before node x can be removed the prop- 
erty of interest is averaged over all the pathways that include the edges between 
nodes x and i G The averaging is performed separately for every node i. 

We will use the waiting time as an example of such a property and show that the 
mean first-passage time in the original and transformed graphs is the same. 

The theory is an extension of the results used to perform jumps to second 
neighbours in previous KMC simulations.^ Consider KMC trajectories that arrive 
at node i, which is adjacent to x. We wish to step directly from i to any vertex 
in the set of nodes F that are adjacent to i or x, excluding these two nodes 
themselves. To ensure that the mean first-passage times from sources to sinks 
calculated in G and G' are the same we must define new branching probabilities, 
from i to all 7 G F, along with a new waiting time for escape from z, r/. 
Here, t[ corresponds to the mean waiting time for escape from i to any 7 G F, 
while the modified branching probabilities subsume all the possible recrossings 
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involving node x that could occur in G before a transition to a node in F. Hence 
the new branching probabilities are:^'^ 



This formula can also be applied if either j or P^^^ vanishes. 

It is easy to show that the new branching probabilities are normalised: 



The modified branching probabilities and waiting times could be used in a KMC 
simulation based upon graph G'. Here we continue to use the notation of §1.4, 
where sinks correspond to nodes a E A and sources to nodes m h E and G 
contains all the nodes in Q expect for the A set, as before. Since the modified 
branching probabilities, P^j, in G' subsume the sums over all paths from i to 
7 that involve x we would expect the sink probability, T^a^ = J^^ea^b^i^ ^ 
trajectories starting at b ending at sink a, to be conserved. However, since each 
trajectory exiting from 7 G F acquires a time increment equal to the average 
value, r/, the mean first-passage times to individual sinks, Ta^b = -^^a,b, are not 
conserved in G' (unless there is a single sink). Nevertheless, the overall mean first- 
passage time to A is conserved, i.e. T^^' = . To prove these results formally 
consider the effect of removing node x on trajectories reaching node i G 



m=0 




(27) 



To calculate t[ we use the method of §1.4: 




(28) 
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from a source node b. The sink probability for a particular sink a is 



oo 




(29) 



and similarly for any other node adjacent to x and any other pathway that revisits 
i other than by immediate recrossings of the type i ^ x ^ i. S' is the ensemble 
of all paths for which probabilistic weights cannot be written in the form defined 
by the last term of the above equation. Hence the transformation preserves the 
individual sink probabilities for any source. 

Now consider how removing node x from each trajectory not included in S' 
affects the mean first-passage time, Ta^b, using the approach of §1.4. 



where the tildes indicate that every branching probability Pq,^/3 is replaced by 
Pa^pe^'^P , as above. The first and last terms are unchanged from graph G in this 




(30) 
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construction, but the middle term, 



J C=0 6ea^, ^3^^ 



E-> « J V ^ P'^ ,xPx,i{^'^i ''"x) ~l~ P'^i,i{j'i ~\~ Pi,xPx,i^x ) V ^ « . 
•^siZ^ n - p. p .)2 i^^^' 

is different (unless there is only one sink). However, if we sum over sinks then 
Tlia&A X]52ea^7 ~ ^ ^' '"^^ '^'^^ simplify the sum over 7 as 

EP^,xPx,ii^Ti ~l~ ''"x) ~l~ P^,i{j'i ~\~ Pi,xPx,iTx^ / \ ^ p/ / fQn\ 

(I- p. P 02 = = M'.i'^i- 

The same argument can be applied whenever a trajectory reaches a node adjacent 
to X, so that = T^', as required. 

The above procedure extends the approach of Bortz, Kalos and Lebowitz^ 
to exclude not only the transitions from the current state into itself but also 
transitions involving an adjacent node x. Alternatively, this method could be 
viewed as a coarse-graining of the Markov chain. Such coarse-graining is accept- 
able if the property of interest is the average of the distribution of times rather 
than the distribution of times itself. In our simulations the average is the only 
property of interest. In cases when the distribution itself is sought, the approach 
described here may still be useful and could be the first step in the analysis of 
the distribution of escape times, as it yields the exact average of the distribution. 

The transformation is illustrated in Fig. 3 for the case of a single source. 
Fig. 3 (a) displays the original graph and its parametrisation. During the first it- 
eration of the algorithm node 2 is removed to yield the graph depicted in Fig. 3 (b). 
This change reduces the dimensionality of the original graph, as the new graph 
contains one node and three edges fewer. The result of the second, and final, 
iteration of the algorithm is a graph that contains source and sink nodes only, 
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with the correct transition probabihties and mean waiting time (Fig. 3 (c)). 

We now describe algorithms to implement the above approach and calculate 
mean escape times from complete graphs with multiple sources and sinks. We 
follow the notation of §1.4 and consider a digraph Qn consisting of Nb source 
nodes, Na sink nodes, and Nj intervening nodes. Qn therefore contains the 
subgraph Gni+Nb- 

The result of the transformation of a graph with a single source b and Na sinks 
using Algorithm 2 is the mean escape time 7^ ' and Na pathway probabilities 
V^, C, & A ^ b. Solving a problem with Nb sources is equivalent to solving A'^^ 
single source problems. For example, if there are two sources bi and 62 we first 
solve a problem where only node bi is set to be the source to obtain q'^'^'+'^B ^^^^ 
the pathway sums from bi to every sink node a E A. The same procedure is then 
followed for 62- 

The form of the transition probability matrix P is illustrated below at three 
stages: first for the original graph, then at the stage when all the intervening 
nodes have been removed (line 16 in Algorithm 2), and finally at the end of the 
procedure: 
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(0 





A^B 





I^I 


I ^ B 
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Each matrix is split into blocks that specify the transitions between the nodes of 
a particular type, as labelled. Upon termination, every element in the top right 
block of matrix P is non-zero. 

Algorithm 2 uses the adjacency matrix representation of graph Q^^, for which 
the average of the distribution of mean first-passage times is to be obtained. 
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For efficiency, when constructing the transition probabihty matrix P we order 
the nodes with the sink nodes first and the source nodes last. Algorithm 2 is 
composed of two parts. The first part (lines 1-16) iteratively removes all the 
intermediate nodes from graph Qn io yield a graph that is composed of sink 
nodes and source nodes only. The second part (lines 17-31) disconnects the source 
nodes from each other to produce a graph with A^^ + A''^ nodes and {Na + A'^)^ 
directed edges connecting each source with every sink. Finally, we evaluate the 
mean first-passage time for each source using the transformed graph. 

The computational complexity of lines 1-16 of Algorithm 2 is 0{Nj + NfNs + 
NfMA + NiNl + NiNsNA). The second part of Algorithm 2 (lines 17-31) scales as 
0{Nq + N'^Na). The total complexity for the case of a single source and for the 
case when there are no intermediate nodes is 0{Nf + NfNA) and 0(A^^-|- A^^A^^), 
respectively. The storage requirements of Algorithm 2 scale as O {{Ni + Nb)"^)- 

We have implemented the algorithm in Fortran 95 and timed it for complete 
graphs of different sizes. The results presented in Fig. 4 confirm the overall cubic 
scaling. The program is GPL-licensed^^ and available online. These and other 
benchmarks presented in this work were obtained for a single Intel® Pentium® 4 
3.00 GHz 512Kb cache processor running under the Debian GNU/Linux operat- 
ing system. The code was compiled and optimised using the Intel® Fortran 
compiler for Linux. 

Finally, we note that the GT method described above is a more general version 
of the method we introduced in our previous publication.™ It can be used to treat 
problems with multiple sources, which may be interconnected, thereby allowing 
us to calculate KMC rate constants exactly. Furthermore, as described in §4, it 
can easily be extended to treat problems where sets of sources and sinks overlap, 
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allowing us to calculate the SS rate constants in equation (9). 

3 Applications to sparse random graphs 

Algorithm 2 could easily be adapted to use adjacency- lists-based data struc- 
tures,^^ resulting in a faster execution and lower storage requirements for sparse 
graphs. We have implemented''^ a sparse-optimised version of Algorithm 2 be- 
cause the graph representations of the Markov chains of interest in the present 
work are sparse. ^° 

The algorithm for detaching a single intervening node from an arbitrary graph 
stored in a sparse-optimised format is given in Algorithm 3. Having chosen the 
node to be removed, 7, all the neighbours [3 G A(i7[7] ^-re analysed in turn, as 
follows. Lines 3-9 of Algorithm 3 find node 7 in the adjacency list of node j3. 
If 13 is not a sink, lines 11-34 are executed to modify the adjacency list of node 
/3: lines 13-14 delete node 7 from the adjacency list of /3, while lines 15-30 make 
all the neighbours a G A(ij[7] of node 7 the neighbours of (3. The sym- 
bol denotes the union minus the intersection of two sets, otherwise known as 
the symmetric difference. If the edge /5 — > a already existed only the branch- 
ing probability is changed (line 21). Otherwise, a new edge is created and the 
adjacency and branching probability lists are modified accordingly (lines 26 and 
27, respectively). Finally, the branching probabilities of node j3 are renormalised 
(lines 31-33) and the waiting time for node (3 is increased (line 34). 

Algorithm 3 is invoked iteratively for every node that is neither a source nor a 
sink to yield a graph that is composed of source nodes and sink nodes only. Then 
the procedure described in §2.1 for disconnection of source nodes (lines 17-31 of 
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Algorithm 2) is applied to obtain the mean escape times for every source node. 
The sparse-optimised version of the second part of Algorithm 2 is straightforward 
and is therefore omitted here for brevity. 

The running time of Algorithm 3 is C)[dc^2i£Adj[c] '^*)' where dk is the degree 
of node k. For the case when all the nodes in a graph have approximately the 
same degree, d, the complexity is 0{d^). Therefore, if there are intermediate 
nodes to be detached, and d is of the same order of magnitude as A^, the cost 
of detaching A^ nodes is (9(A^^). The asymptotic bound is worse than that of 
Algorithm 2 because of the searches through adjacency lists (lines 3-9 and lines 
19-24). If d is sufficiently small the algorithm based on adjacency lists is faster. 

After each invocation of Algorithm 3 the number of nodes is always decreased 
by one. The number of edges, however, can increase or decrease depending on 
the in- and out-degree of the node to be removed and the connectivity of its 
neighbours. If node 7 is not directly connected to any of the sinks, and the 
neighbours of node 7 are not connected to each other directly, the total number 
of edges is increased by d^{3 — d^). Therefore, the number of edges decreases 
(by 2) only when d^ G {1,2}, and the number of edges does not change if the 
degree is 3. For > 3 the number of edges increases by an amount that grows 
quadratically with d^. The actual increase depends on how many connections 
already existed between the neighbours of 7. 

The order in which the intermediate nodes are detached does not change the 
final result and is unimportant if the graph is complete. For sparse graphs, how- 
ever, the order can affect the running time significantly. If the degree distribution 
for successive graphs is sharp with the same average, d, then the order in which 
the nodes are removed does not affect the complexity, which is 0{d^N). If the 
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distributions are broad it is helpful to remove the nodes with smaller degrees 
first. A Fibonacci heap min-priority queue^^'''^ was successfully used to achieve 
this result. The overhead for maintaining a heap is increase-key operations (of 
0{\og{N)) each) per execution of Algorithm 3. Fortran and Python implemen- 
tations of Algorithm 3 algorithm are available online. 

Random graphs provide an ideal testbed for the GT algorithm by providing 
control over the graph density. A random graph, Rn, is obtained by starting 
with a set of nodes and adding edges between them at random. In this work 
we used a random graph model where each edge is chosen independently with 
probability {d) / {N — 1), where {d) is the target value for the average degree. 

The complexity for removal of nodes can then be expressed as 



where dc{i) is the degree of the node, c(i), removed at iteration A(ij[c(i)] is 
its adjacency list, and dj^c{i) is the degree of the jth neighbour of that node at 
iteration i. The computational cost given in Eq. (34) is difficult to express in 
terms of the parameters of the original graph, as the cost of every cycle depends 
on the distribution of degrees, the evolution of which, in turn, is dependent on 
the connectivity of the original graph in a non-trivial manner (see Fig. 5). The 
storage requirements of a sparse-optimised version of GT algorithm scale linearly 
with the number of edges. 

To investigate the dependence of the cost of the GT method on the number of 
nodes, A^, we have tested it on a series of random graphs Rn for different values 
of A^ and fixed average degree, {d). The results for three different values of {d) 
are shown in Fig. 6. The motivation for choosing {d) from the interval [3, 5] was 




(34) 
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the fact that most of our stationary point databases have average connectivities 
for the local minima that fall into this range. 

It can be seen from Fig. 6 that for sparse random graphs Rn the cost scales 
as 0{N'^) with a small ((i)-dependent prefactor. The dependence of the compu- 
tational complexity on (d) is illustrated in Fig. 7. 

From Fig. 5 it is apparent that at some point during the execution of the 
GT algorithm the graph reaches its maximum possible density. Once the graph 
is close to complete it is no longer efficient to employ a sparse-optimised algo- 
rithm. The most efficient approach we have found for sparse graphs is to use the 
sparse-optimised GT algorithm (SGT) until the graph is dense enough, and then 
switch to dense-optimised GT method (DGT), for which pseudocode is given in 
Algorithm 2. We will refer to this approach as SDGT. The change of data struc- 
tures constitutes a negligible fraction of the total execution time. Fig. 8 depicts 
the dependence of the CPU time as a function of the switching parameter Rg. 
Whenever the ratio dc(i)/n{i), where n{i) is the number of nodes on a heap at 
iteration i, is greater than Rg, the partially transformed graph is converted from 
the adjacency list format into adjacency matrix format, and the transformation 
is continued using Algorithm 2. It can be seen from Fig. 5 that for the case of 
a random graphs with a single sink, a single source and 999 intermediate nodes 
the optimal values of Rs lie in the interval [0.07,0.1]. 



4 Overlapping Sets of Sources and Sinks 

We now return to the digraph representation of a Markov chain that corresponds 
to the steady-state results discussed in §1.4. A problem with overlapping sets of 
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sources and sinks can easily be converted into an equivalent problem where there 
is no overlap, and then the GT method discussed in §2.1 and §3 can be applied 
as normal. Hence we can evaluate the SS rate constants in equation (9) in a 
deterministic manner. 

As discussed above, solving a problem with n sources reduces to solving n 
single-source problems. We can therefore explain how to deal with a problem of 
overlapping sets of sinks and sources for a simple example where there is a single 
source-sink i and, optionally, a number of sink nodes. 

First, a new node i' is added to the set of sinks and its adjacency lists are 
initialised to AdjOut[i'] = and Adjln[i'] = Adjln[i]. Then, for every node 
j G AdjOut[i] we update its waiting time as Tj = Tj + Ti and add node j to the set 
of sources with probabilistic weight initialised to Pj^iWi, where Wi is the original 
probabilistic weight of source i (the probability of choosing source i from the set 
of sources). Afterwards, the node i is deleted. 

5 Applications to Lennard-Jones Clusters 
5.1 Oh ^ Ih Isomer isation of LJ38 

We have applied the GT method to study the temperature dependence of the rate 
of interconversion between configurations based on truncated octahedra (Oh) and 
icosahedral (Ih) packing in a 38-atom Lennard-Jones cluster (LJag). The PES 
sample was taken from a previous study^ and contained 1740 minima and 2072 
transition states, not including permutation-inversion isomers. The assignment 
was made in Ref. 5 by solving a master equation numerically to find the eigenvec- 
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tor that corresponds to the smallest magnitude non-zero eigenvalue. As simple 
two-state dynamics are associated with exponential rise and decay of occupation 
probabilities there must exist a time scale on which all the exponential contribu- 
tions to the solution of the master equation decay to zero except for the slowest.^ 
The sign of the components of the eigenvector corresponding to the slowest mode 
was used to classify the minima as It or Oh in character.^ 

The above sample was pruned to remove disconnected minima and create a 
digraph representation that contained 759 nodes with 43 sources, 5 sinks, and 
2639 edges. The minimal, average and maximal degree for this graph were 2, 
3.8 and 84, respectively, and the graph density was 4.6 x 10~^. We have used 
the SDGT algorithm with the switching ratio set to 0.08 to transform this graph 
for several values of temperature. In each of these calculations 622 out of 711 
intermediate nodes were detached using SGT, and the remaining 89 intermediate 
nodes were detached using the GT algorithm optimised for dense graphs (DGT). 

An Arrhenius plot depicting the dependence of the rate constant on temper- 
ature is shown in Fig. 9 (a). The running time for the SDGT algorithm was 
1.8 X 10~^ seconds [this value was obtained by averaging over 10 runs and was 
the same for each SDGT run in Fig. 9 (a)]. For comparison, the timings ob- 
tained using the SGT and DGT algorithms for the same problem were 2.0 x 10^^ 
and 89.0 x 10^^ seconds, respectively. None of the 43 total escape probabilities 
(one for every source) deviated from unity by more than 10"^ for temperatures 
above T = 0.07 (reduced units). For lower temperatures the probability was not 
conserved due to numerical imprecision. 

The data obtained using SDGT method was compared with results from KMC 
simulation, which require increasingly more CPU time as the temperature is 
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lowered. Fig. 9 also shows the data for KMC simulations at temperatures 0.14, 
0.15, 0.16, 0.17 and 0.18. Each KMC simulation consisted of 1000 separate KMC 
trajectories from which the averages were computed. The cost of each KMC 
calculation is proportional to the average trajectory length, which is depicted in 
Fig. 9 (b) as a function of the inverse temperature. The CPU timings for each 
of these calculations were (in order of increasing temperature, averaged over five 
randomly seeded KMC runs of 1000 trajectories each): 125, 40, 18, 12, and 7 
seconds. It can be seen that using the GT method we were able to obtain kinetic 
data for a wider range of temperatures and with less computational expense. 

For the temperatue range of [0.09,0.18] Arrhenius plots for this system were 
obtained previously using KMC and master equation approaches, as well as the 
MM method.^' The results reported here are in quantitative agreement with 
these from Ref. 14. 

5.2 Internal Diffusion in LJ55 

We have applied the graph transformation method to study the centre-to-surface 
atom migration in 55-atom Lennard- Jones cluster. The global potential energy 
minimum for LJ55 is a Mackay icosahedron, which exhibits special stability and 
'magic number' properties. ^^'^^ Centre-to-surface and surface-to-centre rates of 
migration of a tagged atom for this system were considered in previous work.^^'^^ 
In Ref. 14 the standard DPS procedure was applied to create a stationary point 
database based on paths linking the global minimum with the tagged atom occu- 
pying the central position to the same structure with the tagged atom in a surface 
site (There are two inequivalent surface sites lying on two-fold and five-fold axes 
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of rotation). We have reused this database in the present work. 

The sample contained 9907 minima and 19384 transition states. We excluded 
transition states corresponding to degenerate rearrangements from consideration 
because they do not affect the rates. For minima interconnected by more than 
one transition state we added the rate constants to calculate the branching prob- 
abilities. Four digraph representations were created with minimum degrees of 1, 
2, 3 and 4 via iterative removal of the nodes with degrees that did not satisfy 
the requirement. These digraphs will be referred to as digraphs 1, 2, 3 and 4, re- 
spectively. The corresponding parameters are summarised in Table 1. We quote 
CPU timings for the DGT, SGT and SDGT methods for each of these graphs 
in the last three columns of Table 1. Each digraph contained two source nodes 
labelled 1 and 2 and a single sink. The sink node corresponds to the global min- 
imum with the tagged atom in the centre. It is noteworthy that the densities of 
the graphs corresponding to both our examples (LJ38 and LJ55) are significantly 
lower than the values predicted for a complete sample, which makes the use of 
sparse-optimised methods even more advantageous. From Table 1 it is clear that 
the SDGT approach is the fastest, as expected; we will use SDGT for all the rate 
calculations in the rest of this section. 

For this example KMC calculations are unfeasible at temperatures lower than 
about T = 0.3 (reduced units throughout). Already for T = 0.4 the average KMC 
trajectory length is 7.5 x 10^ (value obtained by averaging over 100 trajectories). 
In Ref. 14 it was therefore necessary to use the steady-state approximation for 
the intervening minima to calculate the rate constant at temperatures below 
0.35. These results were not fully converged, and were revised in Ref. 34 using 
alternative formulations of the kinetics. Here we report results that are in direct 
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correspondence with the KMC rate constants, for temperatures as low as 0.1. 

Fig. 10 presents Arrhenius plots for rate constants calculated using the SDGT 
method. The data denoted with circles are the results from seven SDGT calcula- 
tions at temperatures T G {0.3, 0.35, . . . , 0.6} conducted for each of the digraphs. 
The total escape probabilities, and S^, calculated for each of the two sources 
at the end of the calculation deviated from unity by no more than 10~^. For 
higher temperatures and smaller digraphs the deviation was lower, being on the 
order of 10~^° for digraph 4 at T = 0.4. 

At temperatures lower than 0.3 the probability deviated by more than 10~^ 
due to numerical imprecision. This problem was partially caused by the round- 
off errors in evaluation of the terms 1 — -Po,/3-P/3,q, which increase when Pa,pPp,a 
approaches unity. These errors can propagate and amplify as the evaluation 
proceeds. By writing 

Pa,f) = 1 - Pi,f3 = 1 - ea,p 

(35) 

and Pi3^a = 1 — P'y,a = 1 — 6/3^0, 

and then using 

1 — Pa,f3Pl3,a = ^a,p " ^a,()^l3,a + ^P,a (36) 

we were able to decrease 1 — by several orders of magnitude at the expense of 
doubling the computational cost. The SDGT method with probability denomi- 
nators evaluated in this fashion will be referred to as SDGTl. 

Terms of the form 1 — Pa^pPfi^a approach zero when nodes a and (3 become 'ef- 
fectively' (i.e. within available precision) disconnected from the rest of the graph. 
If this condition is encountered in the intermediate stages of the calculation it 
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could also mean that a larger subgraph of the original graph is effectively discon- 
nected. The waiting time for escape if started from a node that belongs to this 
subgraph tends to infinity. If the probability of getting to such a node from any 
of the source nodes is close to zero the final answer may still fit into available 
precision, even though some of the intermediate values cannot. Obtaining the 
final answer in such cases can be problematic as division-by-zero exceptions may 
occur. 

Another way to alleviate numerical problems at low temperatures is to stop 
round-off errors from propagation at early stages by renormalising the branching 
probabilities of affected nodes (3 G Ac(7[7] after node 7 is detached. The corre- 
sponding check that the updated probabilities of node (3 add up to unity could 
be inserted after line 33 of Algorithm 3 (see Appendix 7), and similarly for Al- 
gorithm 2. A version of SDGT method with this modification will be referred to 
as SDGT2. 

Both SDGTl and SDGT2 have similarly scaling overheads relative to the 
SDGT method. We did not find any evidence for the superiority of one scheme 
over another. For example, the SDGT calculation performed for digraph 4 at 
T = 0.2 yielded = TfWi + T^^W2 = 6.4 x 10"^^ and precision was lost 
as both Sf and were less than 10~^. The SDGTl calculation resulted in 

= 8.7 X 10-22 and Sf = = 1.0428, while the SDGT2 calculation produced 

= 8.4 X 10-22 with Sf = = 0.99961. The CPU time required to transform 
this graph using our implementations of the SDGTl and SDGT2 methods was 
0.76 and 0.77 seconds, respectively. 

To calculate the rates at temperatures in the interval [0.1,0.3] reliably we 
used an implementation of the SDGT2 method compiled with quadruple preci- 

32 



sion (SDGT2Q) (note that the architecture is the same as in other benchmarks, 
i.e. with 32 bit wide registers). The points denoted by triangles in Fig. 10 are the 
results from 4 SDGT2Q calculations at temperatures T G {0.10, 0.35, . . . , 0.75}. 
These results are in agreement with those previously reported for this system. 

Using the SDGT2Q formulation we were also able to reach lower temperatures 
for the LJ38 example presented in the previous section. The corresponding data 
is shown in Fig. 9 (triangles). 

6 Conclusions 

The most important result of this work is probably the graph transformation (GT) 
method. The method works with a digraph representation of a Markov chain and 
can be used to calculate the first moment of a distribution of the first-passage 
times, as well as the total transition probabilities for an arbitrary digraph with 
sets of sources and sinks that can overlap. The calculation is performed in a non- 
iterative and non-stochastic manner, and the number of operations is independent 
of the simulation temperature. 

We have presented three implementations of the GT algorithm: sparse-optimised 
(SGT), dense-optimised (DGT), and hybrid (SDGT), which is a combination of 
the first two. The SGT method uses a Fibonacci heap min-priority queue to de- 
termine the order in which the intermediate nodes are detached to achieve slower 
growth of the graph density and, consequently, better performance. SDGT is 
identical to DGT if the graph is dense. For sparse graphs SDGT performs better 
then SGT because it switches to DGT when the density of a graph being trans- 
formed approaches the maximum. We find that for SDGT method performs well 
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both for sparse and dense graphs. The worst case asymptotic scahng of SDGT is 
cubic. 

We have also suggested two versions of the SDGT method that can be used 
in calculations where a greater degree of precision is required. The code that 
implements SGT, DGT, SDGT, SDGTl and SDGT2 methods is available for 
download.''^ 

The connection between the steady-state (SS) kinetic formulation and KMC 
approaches was discussed in terms of digraph representations of Markov chains. 
We showed that rate constants obtained using both the KMC or SS methods can 
be computed using graph transformation. We have presented applications to the 
isomerisation of the L J38 cluster and internal diffusion in the L J55 cluster. Using 
the GT method we were able to calculate rate constants at lower temperatures 
and with less computational expense. 

We also obtained analytic expressions for the total transition probabilities in 
arbitrary digraphs in terms of combinatorial sums over pathway ensembles. It is 
hoped that these results will help in further work, for example, obtaining higher 
moments of the distribution of the first-passage times for arbitrary Markov chains. 

7 Acknowledgements 

S.A.T. is grateful to Cambridge Commonwealth Trust/Cambridge Overseas Trust 
and Darwin College for the financial support. Authors would like to thank Tim 
James, Dwaipayan Chakrabarti and Dr. David J. C. MacKay for comments on 
the manuscript. 



34 



References 

[1] G. Bolch, S. Greiner, H. de Meer and K. S. Trivedi, Queueing Networks and 
Markov Chains, John Wiley and Sons, New York (1998). 

[2] G. R. Grimmett and D. R. Stirzaker, Probability and Random Processes, 
Oxford University Press, Oxford (2005). 

[3] G. R. Grimmett and D. R. Stirzaker, One Thousand Exercises in Probability, 
Oxford University Press, Oxford (2005). 

[4] D.J. Wales, Energy Landscapes: Applications to Clusters, Biomolecules and 
Classes, Cambridge University Press, Cambridge (2003). 

[5] D. J. Wales, Mol. Phys. 100, 3285 (2002). 

[6] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, Elsevier, 
Amsterdam (1981). 

[7] A. B. Bortz, M. H. Kalos and J. L. Lebowitz, J. Comput. Phys. 17, 10 
(1975). 

[8] K. A. Fichthorn and W. H. Weinberg, J. Chem. Phys 95, 1090 (1991). 

[9] M. A. Miller, Energy Landscapes and Dynamics of Model Clusters, Ph.D. 
thesis. University of Cambridge (March 1999). 

[10] M. Block, R. Kunert, E. Scholl, T. Boeck and T. Teubner, New J. Phys. 6, 
166 (2004). 

[11] D. Mukherjee, C. G. Sonwane and M. R. Zachariah, J. Chem. Phys. 119, 
3391 (2003). 

35 



[12] F. M. Bulnes, V. D. Pereyra and J. L. Riccardo, Phys. Rev. E 58, 8692 
(1998). 

[13] R. E. Kunz, Dynamics of First-Order Phase Transitions, Deutsch, Thun 
(1995). 

[14] D. J. Wales, Mol. Phys. 102, 883 (2004). 

[15] D. A. Evans and D. J. Wales, J. Chem. Phys. 121, 1080 (2004). 

[16] D. A. Evans and D. J. Wales, J. Chem. Phys. 119, 9947 (2003). 

[17] M. S. Apaydin, D. L. Brutlag, C. Guestrin, D. Hsu and J.-C. Latombe, in 
RECOMB '02: Proceedings of the Sixth Annual International Conference on 
Computational Biology, pp. 12-21. ACM Press, New York, NY (2002). 

[18] M. S. Apaydin, C. E. Guestrin, C. Varma, D. L. Brutlag and J.-C. Latombe, 
Bioinf. 18, S18 (2002). 

[19] N. Singhal, C. D. Snow and V. S. Pande, J. Chem. Phys. 121, 415 (2004). 

[20] G. Chartrand, Introductory Graph Theory, Dover Publications, New York 
(1977). 

[21] T. H. Gormen, C. E. Leiserson, R. L. Rivest and C. Stein, Introduction to 
Algorithms, The MIT Press, Cambridge, Massachusetts , Cambridge, Mass., 
2nd edn. (2001). 

[22] NIST, Dictionary of Algorithms and Data Structures, 
http:/ /www. nist.gov/dads/ (2005). 

[23] Wikipedia, the Free Encyclopedia, http:/ /www. wikipedia.oig/ (2006). 

36 



[24 
[25; 

[26; 
[27; 



[28' 

[29 
[30 

[31 
[32 
[33 

[34 
[35 
[36 



A. Bar-Haim and J. Klafter, J. Chem. Phys. 109, 5187 (1998). 
M. A. Novotny, Phys. Rev. Lett. 74, 1 (1995). 

M. G. Bulmer, Principles of Statistics, Dover Publications, New York (1979). 

B. E. Trumbo, Relationship Between the Poisson and Exponential Distribu- 
tions, http: / /www. sci.csuhayward.edu/statistics/Resources/Essays/PoisExp. htm 
(1999). 

T. F. Middleton, Energy Landscapes of Model Glasses, Ph.D. thesis. Univer- 
sity of Cambridge (March 2003). 

D. A. Reed and G. Ehrhch, Surf. Sci. 105, 603 (1981). 

A. F. Voter, in Radiation Effects in Solids, pp. 1-22. Springer- Verlag, New 
York (2005). 

M. Raykin, J. Phys. A 26, 449 (1992). 

K. P. N. Murthy and K. W. Kehr, Phys. Rev. A 40, 2082 (1989). 

S. A. Trygubenko, Pathways and Energy Landscapes, Ph.D. thesis. University 
of Cambridge (January 2006). 

S. A. Trygubenko and D. J. Wales, Mol. Phys., in press (2006). 
I. Goldhirsch and Y. Gefen, Phys. Rev. A 33, 2583 (1986). 
I. Goldhirsch and Y. Gefen, Phys. Rev. A 35, 1317 (1987). 



[37] E. W. Weisstein, 'Complete Graph'. From MathWorld — A Wolfram Web 
Resource, http:/ /mathworld. wolfram.com/CompleteGraph. html (2005). 

37 



[38] S. Goldberg, Probability: An Introduction, Dover Publications, New York 
(1960). 



[39 
[40 
[41 
[42 
[43 

[44; 
[45 
[46 

[47; 
[48 
[49 
[50 
[51 
[52 



B. Kahng and S. Redner, J. Phys. A 22, 887 (1989). 

D. Zheng, Y. Liu and Z. D. Wang, J. Phys. A 28, L409 (1995). 

S. K. Kim and H. H. Lee, J. Appl. Phys. 78, 3809 (1995). 

P. C. Bressloff, V. M. Dwyer and M. J. Kearney, J. Phys. A 29, 1881 (1996). 

S. Revathi, V. Balakrishnan, S. Lakshmibala and K. P. N. Murthy, Phys. 
Rev. E 54, 2298 (1996). 

K. Kim, J. S. Choi and Y. S. Kong, J. Phys. Soc. Jpn. 67, 1583 (1998). 

K. Kim, G. H. Kim and Y. S. Kong, Fractals 8, 181 (2000). 

J. Asikainen, J. Heinonen and T. Ala-Nissila, Phys. Rev. E 66, 066706 
(2002). 

P. A. Pury and M. O. Caceres, J. Phys. A 36, 2695 (2003). 

M. Slutsky, M. Kardar and L. A. Mirny, Phys. Rev. E 69, 061903 (2004). 

M. Slutsky and L. A. Mirny, Biophys. J. 87, 4021 (2004). 

D. ben Avraham, S. Redner and Z. Cheng, J. Stat. Phys. 56, 437 (1989). 

Y. Gefen and I. Goldhirsch, Physica D 38, 119 (1989). 

O. Matan and S. Havlin, Phys. Rev. A 40, 6573 (1989). 



[53] H. Haucke, S. Washburn, A. D. Benoit, C. P. Umbach and R. A. Webb, 
Phys. Rev. B 41, 12454 (1990). 

38 



[54] S. Revathi and V. Balakrishnan, J. Phys. A 26, 467 (1993). 

[55] S. H. Noskowicz and I. Goldhirsch, Phys. Rev. A 42, 2047 (1990). 

[56] S. Revathi and V. Balakrishnan, Phys. Rev. E 47, 916 (1993). 

[57] V. Balakrishnan and C. Vandenbroeck, Physica A 217, 1 (1995). 

[58] A. R. Kerstein and R. B. Pandey, Phys. Rev. B 35, 3575 (1987). 

[59] Y. Gefen and I. Goldhirsch, Phys. Rev. B 35, 8639 (1987). 

[60] Y. Gefen and I. Goldhirsch, J. Phys. A 18, L1037 (1985). 

[61] J. W. Haus and K. W. Kehr, Phys. Reports 150, 263 (1987). 

[62] I. G. S. H. Noskowicz, J. Stat. Phys. 48, 255 (1987). 

[63] R. Landauer and M. Buttiker, Phys. Rev. B 36, 6255 (1987). 

[64] R. Tao, J. Phys. A 20, 6151 (1987). 

[65] J. Koplik, S. Redner and D. Wilkinson, Phys. Rev. A 37, 2619 (1988). 

[66] N. M. van Dijk, Queueing Networks and Product Forms, John Wiley and 
Sons, New York (1993). 

[67] E. Gelenbe and G. Pujolle, Introduction to Queueing Networks, John Wiley 
and Sons, New York (1998). 

[68] A. E. Conway and N. D. Georganas, Queueing Networks - Exact Computa- 
tional Algorithms, The MIT Press, Cambridge, Massachusetts (1989). 

[69] D. Eppstein, Z. Galil and G. F. Italiano, in Algorithms and Theory of Com- 
putation Handbook, edited by M. J. Atallah, chap. 8, CRC Press (1999). 

39 



[70] B. V. Cherkassky, A. V. Goldberg and T. Radzik, in SODA '94: Proceedings 
of the Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 
516-525, Society for Industrial and Applied Mathematics, Philadelphia, PA 
(1994). 

[71] G. Ramahngam and T. Reps, J. Algorithms 21, 267 (1996). 



[72 
[73 
[74; 

[75 
[76 
[77 
[78 
[79 

[80 
[81 



G. Ramalingam and T. Reps, Theor. Comput. Sci. 158, 233 (1996). 

GNU General Public License, http:/ /www. gnu.org/copyleft/gpl. html (2006). 

S. A. Trygubenko, Graph Transformation Program, 

http:/ /www. trygub.com/gt/ (2006). 

Debian — The Universal Operating System, http:/ /www. debian.org/ (2006). 

S. A. Trygubenko and D. J. Wales, J. Chem. Phys. 120, 2082 (2004). 

M. L. Fredman and R. E. Tarjan, Journal of the ACM 34, 596 (1987). 

P. Labastie and R. L. Whetten, Phys. Rev. Lett. 65, 1567 (1990). 

R. S. Berry, T. L. Beck, H. L. Davis and J. Jellinek, Adv. Chem. Phys. 70B, 
75 (1988). 

J. P. K. Doye and C. P. Massen, J. Chem. Phys. 122, 084105 (2005). 
J. Pitman, Probability, Springer- Verlag, New York (1993). 



40 



Figure Captions 

1. BKL algorithm for propagating a KMC trajectory applied to a three-state 
Markov chain, (a) The transition state diagram is shown where states and 
transitions are represented by circles and directed arrows, respectively. The 
Markov chain is parametrised by transition probabilities Pa,i, -P/3,1 and Pi 1. 
Absorbing states a and /5 are shaded. If Pi^i is close to unity the KMC 
trajectory is likely to revisit state 1 many times before going to a or j3. (b) 
State 1 is replaced with state 1'. The new Markov chain is parametrised 
by transition probabilities Pa,i', -P/3,1' and the average time for escape from 
state 1 is Ti. Transitions from state 1' to itself are forbidden. Every time 
state 1' is visited the simulation 'clock' is incremented by Ti. 

2. The computational cost of the kinetic Monte Carlo and matrix multiplica- 
tion methods as a function of escape probabilities for K^. M is the number 
of matrix multiplications required to converge the value of the total prob- 
ability of getting from node 1 to nodes 1, 2 and 3: the calculation was 
terminated when the change in the total probability between iterations was 
less than 10^^. The number of matrix multiplications M and the average 
trajectory length (/) can be used as a measure of the computational cost 
of the matrix multiplication and kinetic Monte Carlo approaches, respec- 
tively. The computational requirements of the exact summation method 
are independent of S. Note the logio scale on the vertical axis. 

3. The graph transformation algorithm of §2.1 at work, (a) A digraph with 6 
nodes and 9 edges. The source node is node 1 (white), the sinks are nodes 
a, (3 and 7 (shaded), and the intermediate nodes are 2 and 3 (black). The 
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waiting times and transition probabilities that parametrise the graph are 
given below the diagram, (b) Node 2 and all its incoming and outgoing 
edges are deleted from the graph depicted in (a). Two edges (3^1 and 
/3 <— 3 are added. The parameters for this new graph are denoted by primes 
and expressed in terms of the parameters for the original graph, (c) Node 
3 is now disconnected as well. The resulting graph is composed of source 
and sink nodes only. The total probability and waiting times coincide with 
those of -ftTs, as expected. The new parameters are denoted by a double 
prime. 

4. CPU time needed to transform a dense graph G2N using Algorithm 2 as 
a function of A^. The graph G2N is composed of a subgraph and 
sink nodes. The data is shown for six different cases, including a single 
source, and for sources comprising 20, 40, 60, 90, and 100 percent of the 
nodes in K^-, as labelled. The data for the cases 1 and A^ was fitted as 
5.1 X 10"^ A^'^ and 1.5 x 10~^A^^, respectively (curves not shown). For case 
1 only DetachNode operations were performed while for A^ only Disconnect 
operations were used. Note the logio scale on both axes. 

5. Evolution of the distribution of degrees for random graphs of different ex- 
pected degree, {d) = 5, 10, 15, as labelled. This is a colour-coded projection 
of the probability mass function,^' P{d), of the distribution of degrees 
as a function of the number of the detached intermediate nodes, n. The 
straight line shows P{d,n) for complete graph -ft'iooo- All four graphs con- 
tain a single source, 999 intermediate nodes and a single sink. The trans- 
formation was performed using a sparse-optimised version of Algorithm 2 
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with a Fibonacci-heap-based min-priority queue. It can be seen that as 
the intermediate nodes are detached the density of the graph that is being 
transformed grows. The expected degree of the initial graph determines 
how soon the maximum density will be reached. 

6. CPU time needed to transform a sparse random graph R2N using the GT 
approach described in §2.1 as a function of the number of intermediate 
nodes, A^. R2N is composed of a single source node, N sink nodes and A^ — 1 
intervening nodes. For each value of the data for three different values 
of the expected degree, {d) = 3, 4, 5, is shown, as labelled. Solid lines are 
analytic fits of the form ciV^ where c = 2.3 x 10"", 7.4 x 10"", 1.5 x 10"^° 
for (d) = 3, 4, 5, respectively. The CPU time is in seconds. Note the logio 
scale on both axes. 

7. CPU time needed to transform a sparse random graph R2N using the GT 
approach as a function of the expected degree, (d). The data is shown for 
three graphs with = 500, 750 and 1000, as labelled. R2N is composed of 
a single source node, A^ sink nodes and A^ — 1 intermediate nodes. 

8. CPU time as a function of switching ratio Rg for random graphs of different 
expected degree, {d) = 5, 10, 15, as labelled. All three graphs contain a 
single source, 999 intermediate nodes and a single sink. The transformation 
was performed using the sparse-optimised version of Algorithm 2 until the 
the ratio dc(i)/n{i) became greater than R^. Then a partially transformed 
graph was converted into adjacency matrix format and transformation was 
continued with Algorithm 2. The optimal value of Rg lies in the interval 
[0.07,0.1]. Note the logio scale on both axes. 
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9. (a) Arrhenius plots for the LJ38 cluster, k is the rate constant correspond- 
ing to transitions from icosahedral-like to octahedral- like regions. Green 
circles represent the data obtained from 23 SDGT runs at temperatures 
T G {0.07, 0.075, 0.18}. The data from five KMC runs is also shown 
(squares). The data shown using triangles corresponds to temperatures 
T e {0.035,0.04, ...,0.065} and was obtained using the SDGT2 method 
(discussed in §5.2) with quadruple precision enabled. In all SDGT runs 
the total escape probabilities calculated for every source at the end of the 
calculation deviated from unity by no more then 10"^. For this stationary 
point database the lowest temperature for which data was reported in pre- 
vious works was T = 0.03. (b) The average KMC trajectory length [data in 
direct correspondence with KMC results shown in (a)]. A solid line is used 
to connect the data points to guide the eye. 

10. Arrhenius plots for four digraphs of varying sizes (see Table 1) created 
from a stationary point database for the LJ55 cluster, k is the rate constant 
corresponding to surface-to-centre migration of a tagged atom. Calculations 
were conducted at T G {0.3, 0.35, . . . , 0.7} using the SDGT method (circles) 
and T G {0.1,0.15, . . . ,0.25} using SDGT2Q (triangles). For each of the 
digraphs the calculations yielded essentially identical results, so data points 
for only one of them are shown. Lines connecting the data points are shown 
for each of the digraphs. 
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Tables 



Table 1: Properties of four digraphs corresponding to the LJ55 PES sample from an 
internal diffusion study. \V\ is the number of nodes; \E\ is the number of directed edges; 
d-mim {d) and dmax are the minimum, average and maximum degrees, respectively; p is 
the graph density, defined as a ratio of the number of edges to the maximum possible 
number of edges; r and d are the graph radius and diameter, defined as the maximum 
and minimum node eccentricity, respectively, where the eccentricity of a node v is 
defined as the maximum distance between v and any other node. (/) is the average 
distance between nodes. The CPU time, t, necessary to transform each graph using 
the DGT, SGT and SDGT methods is given in seconds for a sing le 32-bit Intel® 
Pentium® 4 3.00GHz 512Kb cache processor. 
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Algorithms 



Algorithm 1 Calculate the pathway sum S^^ for an arbitrary digraph Gn 

Require: 1 < and a, 6 G {0, 1, 2, . . . , — 1}. ly is a boolean array of size N 
with every element initially set to True. is the number of True elements 
in array W (initialised to A^). P[i,j] is the probability of branching from 
node j to node i. Adjln[i] and AdjOut[i] are the lists of indices of all nodes 
connected to node i via incoming and outgoing edges, respectively. 
Recursive function F{a, (3, W, Nw) 

1: W[/3] ^ False; Nw ^ Nw - I 

2: if a = P and Nw = then 

3: S ^ 1 

4: else 

5: S^O.O 

6: for all i E AdjOut[(3] do 

7: for all j G Adjln[(5] do 

8: if W[i] and W\j] then 

9: E ^ S + P[(3, j]F{j, z, W, Nw)P[i, P] 

10: end if 

11: end for 

12: end for 

13: S^l/(1-S) 

14: if a l3 then 

15: A ^ 0.0 

16: for all i E AdjOut[l3] do 

17: if W[i] then 

18: A + F{a,i,W,Nw)P{hP) 

19: end if 

20: end for 

21: S ^ EA 

22: end if 
23: end if 

24: W[(3] ^ True; Nw ^ Nw + I 
25: return S 
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Algorithm 2 Calculate the total transition probabilities ^ and the mean 
escape times T^'^ in a dense graph 

Require: Nodes are numbered 0, 1, 2, . . . , — 1. Sink nodes are indexed first, 
source nodes last, i is the index of the first intermediate node, s is the index of 
the first source node. If there are no intermediate nodes then i = s, otherwise 
i < s. 1 < N and i,s G {0, 1, 2, . . . , — 1}. T[a] is the waiting time for node 
a, a & {i,i + l,i + 2, . . . , N — 1}. P[i,j] is the probability of branching from 
node j to node i. 



1: for all 7 e {z, -i + 1, z + 2, . . . , s - 1} do 

2: for all /3e {7 + 1,7 + 2,..., A^- 1} do 

3: if P[7,/5] > then 

4: t[P] ^ (r[/3] +r[7]P[7,/3])/(l -P[/3,7]i^[7,/3]) 

5: for all a G {0, 1,2, ...,iV- 1} do 

6: if a 7^ /? and a 7^ 7 then 

7: P[«, (3] ^ {P[a, /3] + P[«, 7]P[7, /3])/(l - P[f3, iMl, P]) 

8: end if 

9: end for 

10: P[7,/3]^0.0 

11: end if 

12: end for 

13: for all a G {0,1,2,..., A^- 1} do 

14: P[a,7]^0.0 

15: end for 

16: end for 

17: for all a G {s, s + 1, s + 2, . . . , A^ - 1} do 

18: for all /5 G {s, s + 1, s + 2, . . . , A^ - 1} do 

19: if a 7^ /3 and P[a, P] > then 

20: P^,f3 ^ P[a,P]; P^,, ^ P[/3,«]; T ^ r[«] 

21: r[«] ^ (rH + r[/3]P^,,)/(l - Pc.,pPp,o.) 

22: t[P] ^ (r[/3] + TP,,;3)/(1 - Pa,pPp,a) 

23: for all 7 G {0, 1, 2, . . . , z - 1} U {s, s + 1, s + 2, . . . , A^ - 1} do 

24: T ^ P[7, a] 

25: P[7, a] ^ (P[7, «] + ^[7, /3]Pl3,a)/il - Pa,f^P(^,a) 

26: P[7, /3] ^ (P[7, /5] + TP,,/3)/(l - 

27: end for 

28: P[a,/3] ^ 0.0; P[/3, a] ^ 0.0 

29: end if 

30: end for 

31: end for 
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Algorithm 3 Detach node 7 from an arbitrary graph Q 



N 



Require: 1 < and 7 G {0, 1, 2, . . . , — 1}. r[z] is the waiting time for node 
i. Adj[i\ is the ordered list of indices of all nodes connected to node i via 
outgoing edges. is the cardinality of v4(ij[z]. 74(i7[z][j] is the index 

of the j'th neighbour of node i. P[i\ is the ordered list of probabilities of 
leaving node i via outgoing edges, \P[i]\ = \Adj[i]\. P[i][j] is the probability 
of branching from node i to node Ac(7[i][j]. 
for all (3^ e {0, 1, 2, ... , \Adj[-f]\ - 1} do 

7/3 ' 1 

for all i G {0, 1,2, . . . , \Adj[l3]\ - 1} do 
if Adj[p][i\ = 7 then 

7/3 ^ i 
break 
end if 
end for 

if not 7^ = — 1 then 

P^,^^l/(l-P[/5][7/3]P[7]N) 

Adj IP] ^ {Adj [P] [0] , Adj IP] [1] , . . . , Adj [P] [7/3 - 1] , Adj [P] [7/3 + 1] , . . . } 

p[p] ^ {p[pm, p[p][i], . . . , p[p][jp - 1], p[/3][7/3 + 1], . . . } 

for all G {0, 1, 2, ... , \Adj['y] \ - 1} do 
a ^ Adj [7] [a^] 
if not a = P then 

if exists edge P —>■ a then 

for all ^ G {0, 1, 2, . . . , \Adj[p]\ - 1} do 
if Ad3\p\\i\ = a then 

P[Pm ^ P[Pm + Pf,,,PMa,] 
break 
end if 
end for 
else 

Adj[P] ^ {a,Adj[pm,Adj[P][llAdj[P][2],...} 
P[P] ^ {P^,,P[7][«,],P[/3][0],P[/3][1],P[ 
end if 
end if 
end for 

for all I G {0,1,2,..., |P[/5]| - 1} do 

p[pm^p[P]m,p 

end for 

T-/3 ^ (T-/3 + Pp,^r^) Pl3,l3 

end if 
end for 
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